function qr_methods_comparison()
% QR_METHODS_COMPARISON QR分解方法比较
% 
% 比较不同QR分解算法的性能和数值特性

fprintf('=== QR分解方法比较 ===\n\n');

%% 1. 算法概述
fprintf('1. QR分解算法概述\n');

algorithms = {
    'Classical Gram-Schmidt', 'CGS', '简单但数值不稳定';
    'Modified Gram-Schmidt', 'MGS', '改进的数值稳定性';
    'Householder反射', 'Householder', '最稳定，广泛使用';
    'Givens旋转', 'Givens', '适合稀疏矩阵';
    'MATLAB内置QR', 'MATLAB', '高度优化的实现'
};

fprintf('QR分解算法对比:\n');
fprintf('算法名称                特点                    适用场景\n');
fprintf('--------------------    --------------------    --------------------\n');
for i = 1:size(algorithms, 1)
    fprintf('%-20s    %-20s    %s\n', algorithms{i, 1}, algorithms{i, 3}, algorithms{i, 2});
end

%% 2. 数值稳定性比较
fprintf('\n2. 数值稳定性比较\n');
fprintf('目标：比较不同算法在病态矩阵上的表现\n\n');

% 测试矩阵：从良态到病态
test_matrices = {
    randn(10, 8),                    % 随机矩阵
    hilb(8),                         % Hilbert矩阵
    create_ill_conditioned(10, 8, 1e12), % 人工病态矩阵
    create_nearly_dependent(12, 8)    % 近似线性相关列
};

matrix_names = {'随机矩阵', 'Hilbert矩阵', '人工病态矩阵', '近似相关矩阵'};

fprintf('数值稳定性测试:\n');
fprintf('矩阵类型        条件数        CGS正交性    MGS正交性    Householder正交性\n');
fprintf('------------    --------      ---------    ---------    -----------------\n');

for i = 1:length(test_matrices)
    A = test_matrices{i};
    cond_A = cond(A);
    
    % Classical Gram-Schmidt
    try
        [Q_cgs, R_cgs] = gram_schmidt_classical(A);
        orthogonality_cgs = norm(Q_cgs' * Q_cgs - eye(size(Q_cgs, 2)), 'fro');
    catch
        orthogonality_cgs = inf;
    end
    
    % Modified Gram-Schmidt
    try
        [Q_mgs, R_mgs] = gram_schmidt_modified(A);
        orthogonality_mgs = norm(Q_mgs' * Q_mgs - eye(size(Q_mgs, 2)), 'fro');
    catch
        orthogonality_mgs = inf;
    end
    
    % Householder QR
    try
        [Q_house, R_house] = householder_qr(A);
        orthogonality_house = norm(Q_house' * Q_house - eye(size(Q_house, 2)), 'fro');
    catch
        orthogonality_house = inf;
    end
    
    fprintf('%-12s    %8.2e    %9.2e    %9.2e    %17.2e\n', ...
            matrix_names{i}, cond_A, orthogonality_cgs, orthogonality_mgs, orthogonality_house);
end

%% 3. 计算复杂度分析
fprintf('\n3. 计算复杂度分析\n');
fprintf('目标：比较不同算法的计算时间\n\n');

sizes = [50, 100, 200, 400];
fprintf('计算时间比较 (秒):\n');
fprintf('矩阵尺寸    CGS        MGS        Householder    Givens     MATLAB\n');
fprintf('--------    -------    -------    -----------    -------    -------\n');

for n = sizes
    m = n;
    k = min(m, n);
    A_timing = randn(m, k);
    
    % Classical Gram-Schmidt
    tic;
    [Q_cgs_time, R_cgs_time] = gram_schmidt_classical(A_timing);
    time_cgs = toc;
    
    % Modified Gram-Schmidt
    tic;
    [Q_mgs_time, R_mgs_time] = gram_schmidt_modified(A_timing);
    time_mgs = toc;
    
    % Householder QR
    tic;
    [Q_house_time, R_house_time] = householder_qr(A_timing);
    time_house = toc;
    
    % Givens QR
    tic;
    [Q_givens_time, R_givens_time] = givens_qr(A_timing);
    time_givens = toc;
    
    % MATLAB内置QR
    tic;
    [Q_matlab_time, R_matlab_time] = qr(A_timing, 0);
    time_matlab = toc;
    
    fprintf('%8d    %7.4f    %7.4f    %11.4f    %7.4f    %7.4f\n', ...
            n, time_cgs, time_mgs, time_house, time_givens, time_matlab);
    
    % 清理变量
    clear Q_*_time R_*_time;
end

%% 4. 内存使用比较
fprintf('\n4. 内存使用比较\n');
fprintf('目标：分析不同算法的内存需求\n\n');

test_size = 500;
A_memory = randn(test_size, test_size);

fprintf('内存使用分析 (矩阵尺寸: %d×%d):\n', test_size, test_size);
fprintf('算法            理论复杂度    实际内存(MB)    相对内存使用\n');
fprintf('------------    ----------    ------------    ------------\n');

% 基准内存使用
base_memory = whos('A_memory');
base_memory_mb = base_memory.bytes / 1024^2;

algorithms_memory = {
    'Classical GS', 'O(mn)', @() gram_schmidt_classical(A_memory);
    'Modified GS', 'O(mn)', @() gram_schmidt_modified(A_memory);
    'Householder', 'O(mn)', @() householder_qr(A_memory);
    'MATLAB QR', 'O(mn)', @() qr(A_memory, 0)
};

for i = 1:size(algorithms_memory, 1)
    alg_name = algorithms_memory{i, 1};
    complexity = algorithms_memory{i, 2};
    alg_func = algorithms_memory{i, 3};
    
    % 测量内存使用
    memory_before = whos;
    total_before = sum([memory_before.bytes]);
    
    [Q_mem, R_mem] = alg_func();
    
    memory_after = whos;
    total_after = sum([memory_after.bytes]);
    
    memory_used = (total_after - total_before) / 1024^2;
    relative_memory = memory_used / base_memory_mb;
    
    fprintf('%-12s    %10s    %12.2f    %12.2f\n', ...
            alg_name, complexity, memory_used, relative_memory);
    
    clear Q_mem R_mem;
end

%% 5. 稀疏矩阵性能
fprintf('\n5. 稀疏矩阵性能\n');
fprintf('目标：比较算法在稀疏矩阵上的表现\n\n');

sparse_sizes = [200, 500];
densities = [0.05, 0.1, 0.2];

fprintf('稀疏矩阵QR分解性能:\n');
fprintf('尺寸    密度    非零元    Householder时间    Givens时间    稀疏优势\n');
fprintf('----    ----    ------    ---------------    ----------    --------\n');

for n = sparse_sizes
    for density = densities
        % 生成稀疏矩阵
        rng(42);
        A_sparse = sprandn(n, round(n*0.8), density);
        A_sparse_full = full(A_sparse);
        nnz_count = nnz(A_sparse);
        
        % Householder QR (稠密)
        tic;
        [Q_house_sparse, R_house_sparse] = householder_qr(A_sparse_full);
        time_house_sparse = toc;
        
        % Givens QR (更适合稀疏)
        tic;
        [Q_givens_sparse, R_givens_sparse] = givens_qr(A_sparse_full);
        time_givens_sparse = toc;
        
        % 稀疏优势
        sparse_advantage = time_house_sparse / time_givens_sparse;
        
        fprintf('%4d    %4.2f    %6d    %15.4f    %10.4f    %8.2f\n', ...
                n, density, nnz_count, time_house_sparse, time_givens_sparse, sparse_advantage);
        
        clear A_sparse A_sparse_full Q_*_sparse R_*_sparse;
    end
end

%% 6. 增量QR分解
fprintf('\n6. 增量QR分解\n');
fprintf('目标：比较批量和增量QR分解\n\n');

base_cols = 100;
additional_cols = [10, 20, 50];
m = 150;

fprintf('增量QR分解性能:\n');
fprintf('基础列数    新增列数    批量QR时间    增量QR时间    加速比\n');
fprintf('--------    --------    ----------    ----------    ------\n');

for add_cols = additional_cols
    % 基础矩阵
    rng(123);
    A_base = randn(m, base_cols);
    [Q_base, R_base] = qr(A_base, 0);
    
    % 新增列
    A_new = randn(m, add_cols);
    A_extended = [A_base, A_new];
    
    % 批量QR分解
    tic;
    [Q_batch, R_batch] = qr(A_extended, 0);
    time_batch = toc;
    
    % 增量QR分解
    tic;
    [Q_incremental, R_incremental] = incremental_qr(Q_base, R_base, A_new);
    time_incremental = toc;
    
    speedup = time_batch / time_incremental;
    
    fprintf('%8d    %8d    %10.4f    %10.4f    %6.2f\n', ...
            base_cols, add_cols, time_batch, time_incremental, speedup);
    
    % 验证结果
    error_incremental = norm(A_extended - Q_incremental * R_incremental, 'fro');
    if error_incremental > 1e-10
        fprintf('警告: 增量QR分解精度损失: %.2e\n', error_incremental);
    end
    
    clear Q_* R_* A_*;
end

%% 7. 并行化潜力
fprintf('\n7. 并行化潜力\n');
fprintf('目标：分析不同算法的并行化特性\n\n');

fprintf('算法并行化特性分析:\n');
fprintf('算法            并行化难度    并行效率    内存访问模式    适合场景\n');
fprintf('------------    ----------    --------    ------------    --------\n');

parallel_analysis = {
    'Classical GS', '困难', '低', '不规则', '教学用途';
    'Modified GS', '困难', '低', '较规则', '中等规模';
    'Householder', '中等', '中等', '规则', '通用';
    'Givens', '容易', '高', '局部化', '稀疏矩阵';
    'Block QR', '容易', '高', '块结构', '大规模'
};

for i = 1:size(parallel_analysis, 1)
    fprintf('%-12s    %10s    %8s    %12s    %s\n', parallel_analysis{i, :});
end

%% 8. 特殊结构矩阵
fprintf('\n8. 特殊结构矩阵\n');
fprintf('目标：测试算法在特殊结构矩阵上的性能\n\n');

% 上Hessenberg矩阵
n_hess = 100;
A_hessenberg = triu(randn(n_hess), -1);  % 上Hessenberg矩阵

fprintf('特殊结构矩阵QR分解:\n');
fprintf('矩阵类型        通用Householder    专用算法    加速比\n');
fprintf('------------    ---------------    --------    ------\n');

% 上Hessenberg矩阵
tic;
[Q_hess_general, R_hess_general] = householder_qr(A_hessenberg);
time_hess_general = toc;

tic;
[Q_hess_special, R_hess_special] = hessenberg_qr(A_hessenberg);
time_hess_special = toc;

speedup_hess = time_hess_general / time_hess_special;

fprintf('%-12s    %15.4f    %8.4f    %6.2f\n', ...
        'Hessenberg', time_hess_general, time_hess_special, speedup_hess);

% 三对角矩阵
A_tridiag = gallery('tridiag', n_hess, -1, 2, -1);

tic;
[Q_tri_general, R_tri_general] = householder_qr(A_tridiag);
time_tri_general = toc;

tic;
[Q_tri_special, R_tri_special] = tridiagonal_qr(A_tridiag);
time_tri_special = toc;

speedup_tri = time_tri_general / time_tri_special;

fprintf('%-12s    %15.4f    %8.4f    %6.2f\n', ...
        '三对角', time_tri_general, time_tri_special, speedup_tri);

%% 9. 数值精度分析
fprintf('\n9. 数值精度分析\n');
fprintf('目标：分析不同算法的数值精度\n\n');

% 使用不同条件数的矩阵
condition_numbers = [1e2, 1e6, 1e10, 1e14];

fprintf('数值精度分析:\n');
fprintf('条件数      MGS重构误差    Householder重构误差    精度比较\n');
fprintf('--------    -----------    -------------------    --------\n');

for target_cond = condition_numbers
    % 构造指定条件数的矩阵
    n_prec = 50;
    [U_prec, ~, V_prec] = svd(randn(n_prec));
    s_prec = logspace(0, -log10(target_cond), n_prec);
    A_prec = U_prec * diag(s_prec) * V_prec';
    
    actual_cond = cond(A_prec);
    
    % Modified Gram-Schmidt
    [Q_mgs_prec, R_mgs_prec] = gram_schmidt_modified(A_prec);
    error_mgs_prec = norm(A_prec - Q_mgs_prec * R_mgs_prec, 'fro') / norm(A_prec, 'fro');
    
    % Householder QR
    [Q_house_prec, R_house_prec] = householder_qr(A_prec);
    error_house_prec = norm(A_prec - Q_house_prec * R_house_prec, 'fro') / norm(A_prec, 'fro');
    
    % 精度比较
    precision_ratio = error_mgs_prec / error_house_prec;
    
    fprintf('%8.0e    %11.2e    %19.2e    %8.2f\n', ...
            actual_cond, error_mgs_prec, error_house_prec, precision_ratio);
    
    clear U_prec V_prec A_prec Q_*_prec R_*_prec;
end

%% 10. 算法选择指南
fprintf('\n10. 算法选择指南\n');

fprintf('QR分解算法选择建议:\n\n');

fprintf('1. 一般用途:\n');
fprintf('   • 推荐: Householder QR (MATLAB内置)\n');
fprintf('   • 原因: 数值稳定，性能优秀，广泛验证\n\n');

fprintf('2. 稀疏矩阵:\n');
fprintf('   • 推荐: Givens旋转\n');
fprintf('   • 原因: 保持稀疏性，局部化操作\n\n');

fprintf('3. 增量更新:\n');
fprintf('   • 推荐: 增量QR算法\n');
fprintf('   • 原因: 避免重复计算，提高效率\n\n');

fprintf('4. 特殊结构:\n');
fprintf('   • Hessenberg矩阵: 专用Hessenberg QR\n');
fprintf('   • 三对角矩阵: 专用三对角QR\n');
fprintf('   • 带状矩阵: 带状QR算法\n\n');

fprintf('5. 教学演示:\n');
fprintf('   • 推荐: Modified Gram-Schmidt\n');
fprintf('   • 原因: 概念清晰，易于理解\n\n');

fprintf('6. 高精度需求:\n');
fprintf('   • 推荐: Householder + 迭代改进\n');
fprintf('   • 原因: 最佳数值稳定性\n\n');

fprintf('7. 大规模并行:\n');
fprintf('   • 推荐: 块QR分解\n');
fprintf('   • 原因: 良好的并行性和缓存效率\n\n');

%% 总结
fprintf('=== QR方法比较总结 ===\n');
fprintf('• Householder QR: 最佳综合性能，推荐用于一般用途\n');
fprintf('• Modified Gram-Schmidt: 概念清晰，适合教学\n');
fprintf('• Givens旋转: 适合稀疏矩阵和并行计算\n');
fprintf('• 专用算法: 针对特殊结构矩阵有显著优势\n');
fprintf('• 算法选择应根据具体应用场景和性能要求\n\n');

end

% 辅助函数实现
function [Q, R] = gram_schmidt_classical(A)
% Classical Gram-Schmidt算法
[m, n] = size(A);
Q = zeros(m, n);
R = zeros(n, n);

for j = 1:n
    v = A(:, j);
    for i = 1:j-1
        R(i, j) = Q(:, i)' * A(:, j);
        v = v - R(i, j) * Q(:, i);
    end
    R(j, j) = norm(v);
    if R(j, j) > eps
        Q(:, j) = v / R(j, j);
    else
        error('矩阵线性相关');
    end
end
end

function [Q, R] = gram_schmidt_modified(A)
% Modified Gram-Schmidt算法
[m, n] = size(A);
Q = A;
R = zeros(n, n);

for j = 1:n
    R(j, j) = norm(Q(:, j));
    if R(j, j) > eps
        Q(:, j) = Q(:, j) / R(j, j);
    else
        error('矩阵线性相关');
    end
    
    for k = j+1:n
        R(j, k) = Q(:, j)' * Q(:, k);
        Q(:, k) = Q(:, k) - R(j, k) * Q(:, j);
    end
end
end

function [Q, R] = householder_qr(A)
% Householder QR分解
[m, n] = size(A);
Q = eye(m);
R = A;

for k = 1:min(m-1, n)
    x = R(k:m, k);
    e1 = zeros(length(x), 1);
    e1(1) = 1;
    
    if norm(x) > eps
        v = sign(x(1)) * norm(x) * e1 + x;
        v = v / norm(v);
        
        % 应用Householder反射
        H_k = eye(m);
        H_k(k:m, k:m) = eye(length(v)) - 2 * (v * v');
        
        R = H_k * R;
        Q = Q * H_k';
    end
end

Q = Q(:, 1:n);
R = R(1:n, :);
end

function [Q, R] = givens_qr(A)
% Givens旋转QR分解
[m, n] = size(A);
Q = eye(m);
R = A;

for j = 1:n
    for i = m:-1:j+1
        if abs(R(i, j)) > eps
            % 计算Givens旋转参数
            a = R(i-1, j);
            b = R(i, j);
            r = sqrt(a^2 + b^2);
            c = a / r;
            s = -b / r;
            
            % 构造Givens矩阵
            G = eye(m);
            G([i-1, i], [i-1, i]) = [c, -s; s, c];
            
            % 应用旋转
            R = G * R;
            Q = Q * G';
        end
    end
end

Q = Q(:, 1:n);
R = R(1:n, :);
end

function A = create_ill_conditioned(m, n, target_cond)
% 创建指定条件数的病态矩阵
[U, ~, V] = svd(randn(m, n));
s = logspace(0, -log10(target_cond), n);
A = U(:, 1:n) * diag(s) * V';
end

function A = create_nearly_dependent(m, n)
% 创建近似线性相关的矩阵
A = randn(m, n);
% 使最后一列近似等于前面列的线性组合
coeffs = randn(n-1, 1);
A(:, end) = A(:, 1:end-1) * coeffs + 1e-12 * randn(m, 1);
end

function [Q, R] = incremental_qr(Q_base, R_base, A_new)
% 增量QR分解
[m, k_base] = size(Q_base);
k_new = size(A_new, 2);

% 投影新列到现有正交空间
P_new = Q_base' * A_new;
A_perp = A_new - Q_base * P_new;

% 对垂直部分进行QR分解
[Q_perp, R_perp] = qr(A_perp, 0);

% 组合结果
Q = [Q_base, Q_perp];
R = [R_base, P_new; zeros(size(R_perp, 1), k_base), R_perp];
end

function [Q, R] = hessenberg_qr(A)
% 上Hessenberg矩阵的专用QR分解
[m, n] = size(A);
Q = eye(m);
R = A;

for k = 1:min(m-1, n)
    if k < m && abs(R(k+1, k)) > eps
        % 只需要一个Givens旋转
        a = R(k, k);
        b = R(k+1, k);
        r = sqrt(a^2 + b^2);
        c = a / r;
        s = -b / r;
        
        % 应用Givens旋转
        G = eye(m);
        G([k, k+1], [k, k+1]) = [c, -s; s, c];
        
        R = G * R;
        Q = Q * G';
    end
end

Q = Q(:, 1:n);
R = R(1:n, :);
end

function [Q, R] = tridiagonal_qr(A)
% 三对角矩阵的专用QR分解（简化实现）
[Q, R] = hessenberg_qr(A);  % 三对角是特殊的Hessenberg
end